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Abstract 

We use the time-dependent variational principle of Balian and Veneroni to derive 
a set of equations governing the dynamics of a trapped Bose gas at finite temper- 
ature. We show that this dynamics generalizes the Gross-Pitaevskii equations in 
that it introduces a consistent dynamical coupling between the evolution of the 
condensate density, the thermal cloud and the "anomalous" density. 
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1 Introduction 



Bose-Einstein condensation (BEC) phenomenon has a long story since its discov- 
ery. Such a fascinating behavior of matter has intrigued many researchers [1]. Static as 
well as dynamic properties of BEC are now intensively studied, both experimentally and 
theoretically, in order to apprehend the way the BEC forms, develops and vanishes. 

Various theoretical techniques have been used successfully, predicting correctly a 
great number of experimental results. Among these techniques, some are most widely 
used, such as the Bogoliubov[2], the Beliaev[3, 4] and the Hartree-Fock-Bogoliubov[5, 6, 7] 
approximations. Some other methods, such as the Monte-Carlo simulations [8] are rising 
interest since they tend to solve the exact quantum many-body problem. 

Although being extremely precise in the static situations, the approximations men- 
tioned above all adopt ad-hoc assumptions about the various quantities such as the order 
parameter <3> (or the condensate density n c ), the non-condensed density or thermal cloud 
n or the anomalous density fh. These assumptions lead inevitably to the fact that the 
approximations are no longer totally consistent in a sense which will become clearer soon. 

In this paper, we rely on a different approach, based on the time-dependent varia- 
tional principle of Balian and Veneroni (BV) [9, 10]. Not only does this method retain the 
essential features of the physics behind the previous approximations, but it also allows 
one to bypass some (if not all) of the ad-hoc assumptions. Indeed, being variational, the 
formalism generates a consistent dynamics for the quantities $, n c , n and m by choosing 
a certain class of trial spaces. 

This well-known advantage of this (and any) variational principle faces however a 
major difficulty related to the choice of the trial spaces. A (difficult) compromise must 
be found between a correct description of the physics on one hand, and the tractability of 
the calculations on the other. In this order of ideas, the BV variational principle requires 
to specify two objects: a density-like operator and a "measured" observable. 

Our choice for the variational spaces consists of a gaussian time-dependent density- 
like operator and a single-particle operator for the observable. This last quantity turns 
out to be of no interest in our particular case but it may play a major role in other 
situations especially when one intends to calculate correlation functions[ll, 12]. 

The machinery we are discussing has in fact already been used elsewhere [10], where 
we have recognized that the variational dynamical equations derived there are a gen- 
eralization of the Gross-Pitaevskii equations[13], that takes into account the coupling 
between the order parameter and the various densities. We called this approach the 
"Time-Dependent Hartree-Fock-Bogoliubov" (TDHFB) approximation. The point is that 
the usual image of a collection of condensed atoms in a bath of thermal particles is not 
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totally true, since the bath has its own dynamics which is sensible to the condensate 
dynamics. 

The purpose of the present paper is to go into further details in this variational 
approximation so as to make a closer connection with other well-known methods (such 
as the ones quoted above) used in the study of Bose-Einstein condensates. Among the 
problems that we intend to study, we can cite in particular the analysis of the static 
properties of the condensate (compared to what is known) as well as the small oscillations 
around the static solutions, when the effects of the coupling with the thermal cloud and 
the anomalous density are taken into account. 

The important paper by A. K. Kerman and P. Tommasini [14] is closely related to 
ours. It uses however the Dirac variational principle which constrains the calculations 
to T = (even if the authors give at the end of the paper a finite temperature pre- 
scription.) The authors also limit themselves to the uniform (that is non trapped) case. 
Therefore, according to what has been said above, we may consider our paper as a twofold 
generalization of [14]. 

The full TDHFB equations also deserve to be solved in order to study the large 
amplitude motion of the condensate and the thermal cloud. Despite its importance, we 
will postpone this study to a future publication which is in progress. 

The paper is organized as follows. In section 2, we recall the major steps used 
in [10, 12] to derive the TDHFB equations and to write them down in the case of the 
BEC problem. Section 3 is devoted to the study of the static solutions, where we recover 
the results of [4] and generalize them to finite temperature. In section 4, we analyze 
the excitations of the condensate by using the RPA technique. Finally, we present some 
concluding remarks and perspectives in order to deal with more complicated situations. 

2 The TDHFB Equations 



The General TDHFB equations were derived in [10] using the BV variational prin- 
ciple. For technical reasons, they were written in terms of the creation and annihilation 
operators a + and a. They may however easily be translated in a more appropriate nota- 
tion for the BEC problem using the boson field operator ^(r) in the Schrodinger picture, 
in exactly the same way as our previous work[12]. Let us first recall some of our notations. 
We introduced the gaussian density operator V(t) (with variational parameters Af(t), \(t) 
and S(t)): 

V(t) = Af(t) exp (A(t)ra) exp (^arS(t)a), (2.1) 
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where a = (a + , a) and r is a symplectic constant matrix. Then, we defined the one-boson 
expectation value (a) and the single-particle density matrix p by: 



(a) 




l + 2p = 



{aa + + a + a) 
2(a+a+) 



—2(aa) 
-{aa + + a + a) 



(2.2) 



which are directly related to X(t) and S(t) respectively (see Eq. A. 6 of ref.[12]). Operators 
such as a are simply the centered operators a — (a). The expectation values are to be 
taken with respect to the density operator (2.1). 

Introducing (2.1) in the BV variational action-like leads, beside the conservation of 
the partition function Z = TrP(t), to what we called the TDHFB equations: 



ih 



d(a) 



= T 



1% 1 



d(H) 
-2 



, d(H) 



(2.3) 



in which (H) is the energy. Some interesting properties are discussed in ([12, 10]). 



In order to make connection with the BEC phenomenon, we introduce first the 
Hamiltonian for trapped bosons [4]: 



H 



h 2 

-— A + V ext (r) -fi 
2m 



a(r) + |^a + (r)a + (r)a(r)a(r), (2.4) 



where 14 xt (r) is the trapping potential, p is the chemical potential and g is the coupling 
constant. The energy S = (H) is easily computed yielding: 

S = J[ -C $ * A$ - * + (Kxt -» + 2gn)\$\ 2 + f |$| 4 - j|nA 

+ (V cxt - p)h + gh 2 + |(|m| 2 + m*$ 2 + m$* 2 )] 



(2.5) 



where the condensate density n c = |$| 2 , the non-condensed density n and the anomalous 
density m are identified respectively with |(a)| 2 , (d + a) and (ad). This is simply because 
we identify in our formalism the boson field operator ^ with the destruction operator a 
and its fluctuation \P with a. 



With these identifications, the Eqs.(2.3) now take the form 

ihh 
ihm 



2 (m*^* -rhM- 
om am 



(2.6) 



which constitutes a closed self-consistent system. The coupling between the order param- 
eter $, the non-condensed density and the anomalous density occurs via the derivatives 
of £ which still contain n and m. 
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Beside the conservation of the energy, the equations (2.6) exhibit the unitary evo- 
lution of the density matrix (already visible in (2.3)) by means of the conservation of the 



We recall the reader that / = coth 2 (hu / kT) for a thermal distribution. 

The expression (2.5) for the energy allows us to write down the Eqs.(2.6) more 
explicitly. They indeed read 



The consistency of our derivation mentioned in the introduction is now clear. Indeed, we 
obtain in Eqs.(2.8) a self-consistent dynamics of the order parameter, the thermal cloud 
and the anomalous density. The equation governing the evolution of $ has been obtained 
elsewhere [4, 15, 16, 17, 18] as an extension of the Gross-Pitaevskii equation, but to our 
knowledge, the two last equations in (2.8), governing the evolution of n and m, were 
never written down before at finite temperature. The exception is ref. [14], where a zero 
temperature (that is / = 1) and uniform (V ext = 0) version was derived. Indeed, taking 
these limits in (2.7) and (2.8) provides the equations obtained in ref. [14]. Therefore, our 
equations are more general. They describe not only a self-consistent dynamics of the non- 
condensed and anomalous densities but also a feedback effect on the order parameter and 
therefore on the condensate density. The coupling is however intimately related to the 
two-body interactions and completely disappears for noninteracting systems, justifying 
therefore the Popov and the Beliaev approximations for weakly interacting atomic gazes. 

It is worth noticing that this dynamics is also number conserving since the total 
density n = n c + n is preserved during the evolution. 

As a final remark, we may note that the information on the temperature is encoded 
in the equation (2.7) which is a property of the density operator (2.1). Indeed, if T is 
specified, then the Heisenberg parameter / (which we recall is a conserved quantity) is cal- 
culable and this in turn allows the computation of n (respectively rh) once m (respectively 
n) is known. The last two equations in (2.8) are therefore not totally independent. 



3 The Static Solutions 

The static properties of the BEC at T = are well known [4]. Indeed, at zero 
temperature, all the atoms are condensed. Therefore, nO = and mO = and n c equals 
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Heisenberg parameter" I defined as p(p + 1) = ^ - 

I = (2n + l) 2 -4|m 



or equivalently 



(2.7) 




(2.8) 



the total density of the gas n. The last two equations in (2.8) become meaningless, and 
the first one gives: 

(-^A + Kxt -fi + gn c ^j $0 = 0. (3.1) 
This provides the density profile in the Thomas- Fermi (TF) approximation [19]: 

n c 0(r) = -(»- Kxt(r)) . (3.2) 
9 

Upon defining the "size" of the fundamental state am = (h/muO) 1/2 and the s-wave 
scattering length a = mg/AixTi 2 , we obtain for a spherical potential V ext (r) = ^muj0 2 r 2 , 
the condensate radius -Rtf and the chemical potential // for a gas of NO bosons as 

R TF =a m (l5N0 7 ^) 1/ \ 
» =^0(l5iV0 z 4) 25 , 

which are well-known expressions in the literature[4, 7]. 

At T > Tbecj there is no condensate so that n c = and nO = n. The static 
TDHFB equations reduce to a single equation which becomes in the TF approximation 

(V ext (r)- f ji + 10gn + g)rM) = 0, (3.4) 

and leads to either a uniform density profile (for rhO = 0), which seems unreasonable in 
the presence of a confining field, or to a non-uniform density of the form 

giving rise to the anomalous density 

W= |„_^zl,| n+ ^,. (3,) 




The maximum spreading of the trapped gas R as well as the chemical potential \i can be 
computed from (3.5). We obtain the results 

R =a HO (l50N0 7 fa) 1/5 , 

/z =s+1^0(l50A^) 2/5 , 

where we can notice the similarities with the totally condensed phase. In the two cases, 
the spreading of the condensate (or of the global system for T > T BEC ) depends essentially 
on the balance between the self-interactions and the confining potential. Furthermore, we 
note on (3.6) that at a maximum distance i? max from the center of the trap 

i? max = iJl-^fc(Vj-l), (3.8) 
V mwO R 
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the anomalous density vanishes. 

When < T < T BEC , we have of course n c ^ and nO ^ 0. In the TF approxima- 
tion, we see that the value mO = can no longer be retained. Therefore, the anomalous 
density, although (maybe) small, is an essential ingredient in the resolution of the static 
equations. 

It is important to notice at this stage that the TF approximation is a somewhat 
hazardous hypothesis for the thermal cloud [20]. Indeed, the traditional image of a con- 
densate surrounded by a smooth thermal cloud is a rather simplified picture. Fortunately, 
we see on Eqs.(2.8) that n can be eliminated in favor of the "relevant " variables n c and 
m. We will henceforth mean by TF approximation the neglect of the kinetic terms in the 
equations of n c and in. 

Let us set £ = i (V cxt (r) — jj) and introduce the parametrization mO = |m0| exp (ia) 
and $0 = V n c exp (i(f>). We then obtain the implicit solutions: 



|m0| 



1 
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nO = 



I 

T=Yq~ 

+ 



;i - 2?) 



"2^3 



n c 



T±Tq + (! - 2«) 
T^ + 3 ( 1 - 2 ^) 



(3.9) 



with g = £ + n. This is obtained by setting a = 20 + 7r which is compatible with the 
Hugenholtz-Pines theorem [21] expressed in our context by the identity |m0| = £ + n + n0. 

The equations (3.9), together with the third static equation in (2.8) may be com- 
bined to yield a quartic equation for q alone which can then be solved numerically to 
provide temperature and position-dependent density profiles. An important preliminary 
result is that the anomalous density is always very small compared to n c or nO whatever 
the conditions are, therefore, justifying a posteriori, the TF approximation used above. 
We shall discuss this and other numerical results in a separate work[22]. 

Nevertheless, one can gain further insights into the static properties at < T < T B ec 
by choosing the parametrization 



1 + 2n0 = Vl cosh a 
2|m0| = vTsinhcr, 



(3.10) 



which automatically satisfies (2.7). Indeed, since we know that the T = case is given 
by a = 0, we may approximately solve the static problem to first order in a. This is 
equivalent to a low temperature expansion (but far from the transition). After some 
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algebra, we obtain 

n c =-S-(y/l-l) + y/lri 

h0 = ^1 + ^^ (3.11) 

|m0| =-^77(1 + 277), 

rj being a small expansion parameter. What we observe on (3.11) is that it is a natural 
extension to finite temperature of an expression like (3.2). To lowest order (77 = 0), the 
result is a temperature-dependent shift of the condensate density with respect to the 
T = case. 

Let us see what happens to the condensate radius and the chemical potential which 
were given by (3.3). We define the condensate radius in the TF approximation by the 
point where n c vanishes. This gives 

V cxt (R TF ) = fi-(Vl-l)g. (3.12) 

The number of condensed atoms N c , which is a measurable quantity, writes 

N c = An J RTF n c 0(r)r 2 dr. (3.13) 

After integrating and using (3.12), we obtain the following remarkable expressions: 

R TF =a HO (l5N c7 ^) 1/ \ 



H =^( v / 7-l) + ^0(l5iV c ^) 



2/5 (3-14) 



The condensate radius is thus given by the same formula as its zero temperature coun- 
terpart, the sole difference lying in the appearance of N c instead of the total number NO. 
For the chemical potential, we observe, as for n c 0, a temperature-dependent shift with 
respect to the T = case, but here also, it is the number of condensed atoms which is 
involved and not the total number of atoms. 

4 Excitations of the Condensate 



The small excitations of the condensate (collective modes) are well studied with the 
RPA technique. Indeed, one may first set 

$ = $0 + 5$ 

n =n0 + 5h (4.1) 
m = mO + 5m, 

$0, nO and mO being the static solutions satisfying (3.1), (3.4) or (3.9) according to the 
phase the system is in and <5$, 5n and 5m are small deviations from local equilibrium. 
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Then, one may expand the dynamical equations (2.8) dropping terms up to second order, 
but keeping the kinetic terms since they are crucial for the excitation spectrum. 

At zero temperature, we recall that the sole meaningful equation is the first one 
in the eqs.(2.8). Its expansion around the static solution (3.1) leads to the following 
equation: 
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2gm 



A + £ + 2n c <5$ + $0 2 5$*. 



(4.2) 



Upon setting e k = h 2 k 2 /2m, we can solve for the modes <5$(r, t) = e lu)t u(r) + e wt v (r) to 
obtain the dispersion relation for a uniform gas 



huj k = Je k (e k + 2gn c 0), 



(4.3) 



which is the Bogoliubov spectrum at T = [7]. For instance, at small momenta, one 
obtains the phonon spectrum u k = \f^^k. 

For finite temperatures and above the transition, we linearize the TDHFB eqs.(2.8) 
without using the TF approximation in order to keep the kinetic term e' k = e k /g. The 
procedure provides a (4 x 4) RPA system of the form: i^V = MV, where V is the column 
vector ((5$, <5$*, 5m, 8m*) and the RPA matrix is given by: 



M = 



( My 

-Ml 
M 4 
V M* 



M 2 
-M x 
-M 5 

—M? 



M* 


M 6 




\ 

-M 3 



-Mj 



(4.4) 



with 



Mi 
M 2 
M 3 
M 4 
M 5 
M 6 



e k + £ + 2n0 
mO - $0 2 
$0 

2$0(1 + 2n0 - n c 0) - 2m0$0* 

2$0(m0 + $0 2 ) 

A(t k + f + 2n) + 1 + 2n0. 



(4.5) 



After some algebra, we obtain the expressions for the eigenfrequencies: 

1/2 



{1 N J 



(4.6) 



where 



and 



B = 4n c 0(l + 2n0) + [4(M X + 2n c 0) 



2n] [4(Mi + n c 0) 



2n], 



E = -8n c [1 + 2n0 + 3M X + 4n c 0] [M x (l + 2n0) + 4(Mi + n c 0)(Mi + 2n c 0)] 
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The static solutions are computed from (3.9) with £ replaced by e' k + £. 

The spectrum (4.6) clearly exhibits a departure from the Bogoliubov spectrum (4.3). 
It is indeed temperature and position-dependent. The two modes uj± reduce respectively 
to efc + V cxt — f-i and 4(efc + V cxt — /i) in the noninteracting case (g — 0). Therefore, one 
may say in the interacting case that and uj + describe the coupled oscillations of the 
condensate fraction and the anomalous density. What remains is a comparison between 
these modes and existing data. But this requires a detailed knowledge of the static solution 
as shown by (4.6). Indeed, the temperature dependence of the static solution complicates 
the expressions for the eigenfrequencies and one therefore cannot use the simplifications 
that appear in the T = case, as was performed in [6, 14]. 

The detailed study of the static solution as well as the eigenmodes is in progress 
and we will postpone the results to a future paper [22]. 

5 Conclusions and Perspectives 

We have been concerned in this paper with a dynamical variational generalization 
of the Gross-Pitaevskii equation, which takes into account the coupling of the condensate 
with the thermal cloud and with the anomalous density. We show that our derivation is 
consistent with all known approximations which go beyond the Gross-Pitaevskii approx- 
imation, namely, the Popov, the Beliaev and the HFB approximations. 

The equations that we obtain are fully self-consistent, mainly because they not only 
introduce a dynamics of the thermal cloud and the anomalous density, but they allow also 
for a consistent feedback effect of these densities on the condensate fraction. 

Instead of solving the full dynamical equations, we choose to focus first on the static 
situation, where theoretical works as well as experimental data do exist. 

At zero temperature, we obtain familiar expressions for the chemical potential and 
the condensate radius. But for finite temperatures and above the transition, the situation 
is much more complicated and requires a numerical study. 

The preliminary results show a good qualitative agreement with what is known; in 
particular, the anomalous density is always extremely small which is compatible with the 
Thomas-Fermi approximation. 

We then turn to the small amplitude motion and derive RPA-like equations which 
provide two coupled modes of oscillations; the well-known breathing modes of the conden- 
sate. A direct comparison with the results of [6, 14] is unfortunately a little bit delicate 
since both the T = and the V ext = cases were considered there. We have shown in 
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particular that, owing to a temperature dependence from the beginning, makes the T = 
limit a rather subtle question. 

This paper is dedicated to the memory of Dominique Vautherin (Arthur) , an active 
member of the Division de Physique Theorique, Institut de Physique Nucleaire (IPN), 
Orsay-France. 

We are grateful to P. Schuck and M. Veneroni for fruitful discussions. We are 
particularly indebted to C. Martin for her valuable comments and a careful reading of the 
manuscript. 
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